## ----------------------------------------------------------------------
##
## Draw maps: Figure 1, 5, A3
##
## ----------------------------------------------------------------------

## environment setting
Sys.setenv(LANGUAGE="en")
gc();gc()
rm(list = ls())
options(scipen = 999)   ## Disable scientific notation

## ----------------------------------------------------------------------
##
## Initial settings
##
## ----------------------------------------------------------------------
## Set directory
## ----------------------------------------------------------------------
root_dir = "YOUR_DIRECTORY/Replication_Folder"
## Load packages and functions
source(file.path(root_dir, "2_Code/1_PackagesFunctions.R"))
needs(lubridate)

## ----------------------------------------------------------------------
## Data directory
data_dir = root_dir	%>%	file.path("1_Data")	%>%	dir_create()

## Output directory
output_dir = root_dir	%>%	file.path("3_Result")	%>%	dir_create()
figure_dir = root_dir	%>%	file.path("4_Figures")	%>%	dir_create()
## Working directory
working_dir = root_dir	%>%	file.path("5_Working")	%>%	dir_create()
## ----------------------------------------------------------------------

## ------------------------------------------------------------
## Load and prepare data objects
## ------------------------------------------------------------
## Census data
tokyo2keep = data_dir	%>%
  file.path("census_data4.csv") %>%
  read_csv(col_types = cols()) %>%
  mutate(group = "FullSample")

## Full polygons
tokyo_full_poly = data_dir %>% file.path("RaidShp_May2020b.rds") %>% read_rds()
background_poly = tokyo_full_poly   %>%
  filter(
    row_id != 941 & row_id != 974 &
      row_id != 1302 & row_id != 168
  )	%>%
  st_union()	%>%
  as("Spatial")
## Raid polygons w/ prewar census population density
raid_poly = tokyo_full_poly %>%
  filter(!is.na(ratio_damage))	%>%
  as("Spatial")
## ------------------------------------------------------------



## ------------------------------------------------------------
##
## Figure 1: Draw residential ratio map
##
## ------------------------------------------------------------
## Target locations
target_locationz = read_rds(file.path(data_dir, "TargetLocationsCombined.rds"))
## ------------------------------------------------------------
## Pull out the imperial palace polygon
palace = tokyo_full_poly    %>%
  filter(row_id == 88)    %>%
  st_centroid %>%
  st_transform(3095)  ## EPSG:3095 --- http://spatialreference.org/ref/epsg/tokyo-utm-zone-54n/
## Buffer of 10km
palace_buffer = st_buffer(palace, dist = 10*10^3)	%>%
  as("Spatial")  %>%
  spTransform(prj)

## Buffer of 10km
palace_buffer_2km = st_buffer(palace, dist = 2*10^3)   %>%
  as("Spatial")   %>%
  spTransform(prj)
## ------------------------------------------------------------
## Set polygon colors
## ------------------------------------------------------------
background_poly_col = "ivory2";
bord_col = "gray10";
residential_colz = c("white","gray40")  ## gray scale
nn = 10
fxBrks = seq(0, 1, by = 0.1)
nn = length(fxBrks)
palz = colorRampPalette(residential_colz, space = "rgb")(nn)
raid_poly$col = palz[raid_poly$ratio_residential*10+1]
circle_col = c("orangered2", "dodgerblue2", "white")
dot_col = "chartreuse2"
  
## ------------------------------------------------------------
## Plot
## ------------------------------------------------------------
fig_height = 6
# dev.new(height = fig_height, width = 8)
png(file.path(figure_dir, "fig_1"), width = 6.15, height = fig_height, units = "in", res = 240)
par(mar=rep(0,4), oma=rep(0,4), bg="white", family = "Helvetica", lend="square")
plot(background_poly, col = background_poly_col, border = bord_col, lwd = .25, xaxs = "i", yaxs = "i") ## Background
plot(raid_poly, col=raid_poly$col, border = bord_col, lwd=.35, add = TRUE)       ## square w/ colors
## 10km Buffer
lines(palace_buffer, lwd = 3, col = circle_col[3])
lines(palace_buffer, lwd = 1, col = circle_col[1])
## 2km Buffer
lines(palace_buffer_2km, lwd = 3, col = circle_col[3])
lines(palace_buffer_2km, lwd = 1, col = circle_col[2])
# text(coordinates(palace_buffer_2km)[,1], y = coordinates(palace_buffer_2km)[,2], labels = "r=2km", col = "maroon2")
## Target locations
points(target_locationz, pch = 21, bg = adjustcolor(dot_col, alpha = 0.85), col = "gray10", cex = 1)
## Add legend
legend(
  title = "Residential Area (%)",
  "bottomright",
  inset=c(.001,.001),
  bg = "white", box.lwd = 1,
  ncol=2, cex=1,
  col = c(palz, background_poly_col),
  legend = c(str_c(c(seq(0, 100, by = 10)), "%"), "NA"),
  lwd=12, seg.len=.85
)
## Close
dev.off()
## ------------------------------------------------------------




## ------------------------------------------------------------
##
## Figure 5a: Draw raid damage map
## Raw human coding
##
## ------------------------------------------------------------
## Set some colors
# background_poly_col = "gray90"
background_poly_col = "ivory2";
sea_col = "azure3";
land_col = "gray96";
bord_col = "gray10";
fig_height = 6
## ------------------------------------------------------------
## Set polygon colors
damage_colz = c("white","orange","firebrick2","darkred")
nn = 10
fxBrks = seq(0, 1, by = 0.1)
nn = length(fxBrks)
palz = colorRampPalette(damage_colz, space = "rgb")(nn)
raid_poly$col = palz[raid_poly$ratio_damage*10+1]
## ------------------------------------------------------------
## Plot base map
## ------------------------------------------------------------
# dev.new(height = fig_height, width = 8)
png(file.path(figure_dir, "fig_5a.png"), width = 6.15, height = fig_height, units = "in", res = 240)
par(mar=rep(0,4), oma=rep(0,4), bg="white", family = "Helvetica", lend="square")
plot(background_poly, col = background_poly_col, border = bord_col, lwd = .25, xaxs = "i", yaxs = "i") ## Background
plot(raid_poly, col=raid_poly$col, border = bord_col, lwd=.35, add = TRUE)       ## square w/ colors
## Scalebar
# scalebar_rev(loc = c(139.585, 35.54), length = .09)
## Add legend
legend(
  title = "Raid Damage (%)",
  "bottomright",
  inset=c(.001,.001),
  bg = "white", box.lwd = 1,
  ncol=2, cex=1,
  col = c(palz, background_poly_col),
  legend = c(str_c(c(seq(0, 100, by = 10)), "%"), "NA"),
  lwd=12, seg.len=.85
)
## Close
dev.off()
## ------------------------------------------------------------



## ----------------------------------------------------------------------
##
## Figure 5b: Draw raid damage map
## Residualized human coding
##
## ----------------------------------------------------------------------
## Prepare data and residualized raid scores
## ----------------------------------------------------------------------
## Treatment
treatment_vec = "ratio_damage"
## ----------------------------------------------------------------------
## X: Pretreatment covariates
## ----------------------------------------------------------------------
## Distance variables
distance_termz_base = c("palace_dist_ln", "minTargetDistance_ln")
## Polynomial
polynomial_degree = 5
distance_termz_poly = polynomial_varname(distance_termz_base, degree = polynomial_degree)
distance_term_specification = list( distance_termz_poly, str_c("s(", distance_termz_base, ")") )
## Geographical features
spatial_covariates = c(
  ## Residential ratio and geographical area
  "ratio_residential", "poly_area_ln", "n_neighbor",
  ## Prewar population density
  "PopDensity_1939_km2_ln",
  ## Terrain variables
  "mean_elevation_ln", "mean_slope_ln", "river_distance_ln",
  ## Railway
  "railway_length_ln", "n_stations", "SpLag_railway_length_ln", "SpLag_n_stations"
)
## Flammability score
flame_scorez = c("hazard", "max_nghb_hazard")
## ----------------------------------------------------------------------
## Spatial spline/polynomials and fixed effects
## ----------------------------------------------------------------------
## District FEs
district_fe = "as.factor(prewar_district)"	## 35-district FE
## Spatial terms
spatial_term_lbl = c("Polynomial", "GAM")   ## This indicates polynomial or spline specification in the csv files
lonlat_base = c("z_lon", "z_lat")
spatial_polynomial = c(
  polynomial_varname(lonlat_base, degree = polynomial_degree),
  "I(z_lon*z_lat)", "I(z_lon^2*z_lat)", "I(z_lon*z_lat^2)","I(z_lon^2*z_lat^2)", "I(z_lon^3*z_lat)", "I(z_lon*z_lat^3)",
  "I(z_lon^4*z_lat)", "I(z_lon*z_lat^4)", "I(z_lon^3*z_lat^2)", "I(z_lon^2*z_lat^3)"
)
spatial_term_specification = list(spatial_polynomial, "te(longitude, latitude)")

## ----------------------------------------------------------------------
## Prepare models
## ----------------------------------------------------------------------
## 1. Polynomial model
remove_polynomial = chr2fml(treatment_vec,
                            idv_list = list(
                              spatial_covariates,
                              distance_term_specification[[1]], spatial_term_specification[[1]],
                              district_fe	## ADDED: Prewar 35-district FEs
                            )
)
## 2. GAM with additional spline terms
remove_gam = chr2fml(treatment_vec,
                     idv_list = list(
                       spatial_covariates,
                       distance_term_specification[[2]], spatial_term_specification[[2]],
                       district_fe	## ADDED: Prewar 35-district FEs
                     )
)
## ----------------------------------------------------------------------
## Pull out variable names
xvarz = fm2xvarWspline(remove_polynomial)
xvarz = xvarz[!str_detect(xvarz, c("I\\(|\\)"))]
varz2keep = c(
  "key_code", "year", treatment_vec,
  # depvar_vec_combined,
  xvarz,
  "palace_dist", "pop", "longitude", "latitude"
)	%>%
  unique
## ----------------------------------------------------------------------
## Extract variables and add subsample index
tokyo2keep = tokyo2keep	%>%
  filter(year == 2010)	%>%
  drop_na(ratio_damage, pop)	%>%
  ## Standardize for polynomial terms
  mutate_at(vars(all_of(distance_termz_base)), list(scale2))

## ----------------------------------------------------------------------
## Note: gam.residuals() function
## --- Response residuals are the raw residuals (data minus fitted values)
## --- Scaled Pearson residuals are raw residuals divided by the standard deviation of the data according to the model mean variance relationship and estimated scale parameter.
## --- Pearson residuals are the same, but multiplied by the square root of the scale parameter
## --- Deviance residuals simply return the deviance residuals defined by the model family.
## --- Working residuals are the residuals returned from model fitting at convergence.
## ----------------------------------------------------------------------
## Estimate models to obtain residuals
## --- WITH 35-district FEs
remove_wls_2010 = lm(remove_polynomial, data = tokyo2keep, weights = tokyo2keep$pop)
remove_gam_2010 = gam(remove_gam, data = tokyo2keep, weights = tokyo2keep$pop)
tokyo2keep$res_ply_raid_wFE = residuals(remove_wls_2010)
tokyo2keep$res_gam_raid_wFE = residuals(remove_gam_2010, type = "response")
## ----------------------------------------------------------------------
tokyo_resid = tokyo2keep	%>%	select(key_code, res_ply_raid_wFE, res_gam_raid_wFE)

# Extract the data frame from the SpatialPolygonsDataFrame
raid_poly_df <- as.data.frame(raid_poly) %>%
  mutate(key_code = as.numeric(as.character(key_code)))

# Perform the join
raid_poly_resid_df <- raid_poly_df %>%
  left_join(tokyo_resid, by = "key_code") 

# Create a new SpatialPolygonsDataFrame with the joined data
raid_poly_resid <- SpatialPolygonsDataFrame(raid_poly, data = raid_poly_resid_df)

## ----------------------------------------------------------------------
## Map
## ----------------------------------------------------------------------
## Set baseline polygon colors
background_poly_col = "ivory2";
sea_col = "azure3";
land_col = "gray92";
bord_col = "gray10";
nn = 11
fig_height = 6
varz_vec = c("res_gam_raid_wFE")
fig_lbl = c("gam_residual_w35dFE")
map_dir = figure_dir
## ----------------------------------------------------------------------
## Residualized raid index
var2plot = raid_poly_resid@data	%>%	pull(varz_vec[1])
var2plot = rescaleIt01(var2plot)	## Rescale to a [0,1] interval
## ----------------------------------------------------------------------
## Color intervals
fxBrks = seq(min(var2plot, na.rm = TRUE), max(var2plot, na.rm = TRUE), length = nn)
palz = colorRampPalette(c("white","orange","firebrick2","darkred"), space = "rgb")(nn)
raid_poly_resid$col = palz[var2plot*10+1]
## ----------------------------------------------------------------------
## Plot map
# dev.new(width = 6.15, height = fig_height)
png(file.path(map_dir, "fig_5b.png"), width = 6.15, height = fig_height, units = "in", res = 240)
par(mar = rep(0,4), oma = rep(0,4), bg="white", family = "Helvetica", lend="square")
# plot(background_poly, col = "gray90", border = NA, lwd = .35)    ## square w/ colors
plot(background_poly, col = background_poly_col, border = bord_col, lwd = .25, xaxs = "i", yaxs = "i")	## Background
plot(raid_poly_resid, col = raid_poly_resid$col, border = bord_col, lwd =.25, add = TRUE)				## Square w/ colors
## Add legend
legend(
  title = "Raid Damage (%)",
  "bottomright",
  inset = c(.001,.001),
  bg = "white", box.lwd = 1,
  ncol=2, cex=1,
  col = c(palz, background_poly_col),
  legend = c(str_c(c(seq(0, 100, by = 10)), "%"), "NA"),
  lwd = 12, seg.len = .85
)
dev.off()
## ----------------------------------------------------------------------



## ------------------------------------------------------------
##
## Figure A3: Prewar population density
##
## ------------------------------------------------------------
## Set some colors
background_poly_col = "ivory2";
sea_col = "azure3";
land_col = "gray96";
bord_col = "gray10"
  ## ------------------------------------------------------------
## Set Colors: Prewar Population density
## ------------------------------------------------------------
fxBrks = c(18, 10^3, 5*10^3, 10^4, 3*10^4)
nn = length(fxBrks)
classes_km = classIntervals(
  raid_poly$PopDensity_1939_km2, style = "fixed", n = nn,
  fixedBreaks = c(fxBrks, max(raid_poly$PopDensity_1939_km2))
)
# classes_km$brks = round(classes_km$brks, 2)
classes_km$brks = round(classes_km$brks)
pal = colorRampPalette(brewer.pal(nn, "BuGn"))(nn)
# pal = palz[[1]]   ## customized palette:
# pal = pal[2:length(pal)]    ## drop white
cols = findColours(classes_km, pal)    ## Set cols
names(attr(cols,"table"))[length(names(attr(cols,"table")))] = str_c("[", fxBrks[length(fxBrks)], ",]")

## ------------------------------------------------------------
## Draw Palace Distance
## ------------------------------------------------------------
png(file.path(figure_dir, "fig_a3.png"), width = 6.15, height = fig_height, units = "in", res = 240)
par(mar=rep(0,4), oma=rep(0,4), bg="white", family = "Helvetica", lend="square")
plot(background_poly, col = background_poly_col, border = bord_col, lwd = .25, xaxs = "i", yaxs = "i") ## Background
plot(raid_poly, col=cols, border = bord_col, lwd=.25, add = TRUE)     ## square w/ colors
## Add legend
legend(
  title = "Population Density (km2, 1939)",
  "bottomright",
  inset=c(.001,.001),
  bg = "white", box.lwd = 1,
  ncol = 1, cex=1,
  col=attr(cols,"palette"),
  legend=names(attr(cols,"table")),
  lwd=12, seg.len=.85)

dev.off()
## ------------------------------------------------------------

##EOF

